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ABSTRACT 

We present numerical simulations and explore scalings and anisotropy of compressible 
magnetohydrodynamic (MHD) turbulence. Our study covers both gas pressure dom- 
inated (high P) and magnetic pressure dominated (low f3) plasmas at different Mach 
numbers. In addition, we present results for superAlfvenic turbulence and discuss in 
what way it is similar to the subAlfvenic turbulence. We describe a technique of sep- 
arating different magnetohydrodynamic (MHD) modes (slow, fast and Alfven) and 
apply it to our simulations. We show that, for both high and low (3 cases, Alfven and 
slow modes reveal the Kolmogorov k"^^^ spectrum and scale-dependent Goldreich- 
Sridhar anisotropy, while fast modes exhibit k~^/'^ spectrum and isotropy. We discuss 
the statistics of density fluctuations arising from MHD turbulence at different regimes. 
Our findings entail numerous astrophysical implications ranging from cosmic ray prop- 
agation to gamma ray bursts and star formation. In particular, we show that the rapid 
decay of turbulence reported by earlier researchers is not related to compressibility and 
mode coupling in MHD turbulence. In addition, we show that magnetic field enhance- 
ments and density enhancements are marginally correlated. Addressing the density 
structure of partially ionized interstellar gas on AU scales, we show that the viscosity- 
damped regime of MHD turbulence that we earlier reported for incompressible flows 
persists for compressible turbulence and therefore may provide an explanation for 
those mysterious structures. 
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C3 1 INTRODUCTION 

Astrophysical turbulence is ubiquitous and it holds the key 
to many astrophysical processes (star formation, heating of 
the interstellar medium, properties of accretion disks, cosmic 
ray transport etc). Therefore understanding of turbulence is 
a necessary requirement for making further progress along 
any of those directions. Unlike laboratory turbulence astro- 
physical turbulence is magnetized and highly compressible. 

Turbulence has been studied in the context of the inter- 
stellar medium (ISM) and the solar wind. The ISM in the 
Milky Way and neighboring galaxies is known to be turbu- 
lent on scales ranging from AUs to kpc (Armstrong, Rick- 
ett & Spangler 1995; Stanimirovic & Lazarian 2001; Desh- 
pande et al. 2000). The solar wind also exhibits small-scale 
turbulence (Leamon et al. 1998). The measured statistics 
of fluctuations in the ISM and the solar wind is consistent 
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with the Kolmogorov turbulence obtained for incompressible 
unmagnetized fluid. This surprising observational evidence^ 
resulted in numerous attempts to use Kolmogorov statistics 
for practical computations of astrophysical quantities, e.g. 
cosmic ray scattering. We shall show below that in most 
cases such sort of calculations brings erroneous results. 

Why would we expect astrophysical fluids to be turbu- 
lent and how can we study astrophysical turbulence? A fluid 
of viscosity v gets turbulent when the rate of viscous dissi- 
pation, which is ~ at the energy injection scale L, is 



^ The ambiguities of the data interpretation were frequently 
quoted to justify ignoring this fact. For instance, electron den- 
sity fluctuations discussed in Armstrong et al. (1995) provide 
only indirect evidence for Kolmogorov-type spectrum. However, 
the solar wind observations are made in situ and more difficult 
to disregard. Moreover, a newly developed statistical technique 
(Lazarian & Pogosyan 2000) allowed us to measure Kolmogorov 
type spectrum of velocity fluctuations (see also Lazarian & Es- 
quivel 2003). 
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much smaller than the energy transfer rate ^ Vl/L, where 
Vl is the velocity dispersion at the scale L. The ratio of the 
two rates is the Reynolds number Re = VlL/v. In general, 
when Re is larger than 10 — 100 the system becomes turbu- 
lent. Chaotic structures develop gradually as Re increases, 
and those with Re ~ 10"^ arc appreciably less chaotic than 
those with Re ~ 10**. Observed features such as star forming 
clouds are very chaotic for Re > 10*. This not only ensures 
that the fluids are turbulent but also makes it difficult to 
simulate the turbulence. The currently available 3D simu- 
lations for 512 cubes can have Re up to 0(10'*) and are 
limited by their grid sizes. Therefore, it is essential to find 
"scaling laws" in order to extrapolate numerical calculations 
(Re ~ O(IO^)) to real astrophysical fluids [Re > 10**). Wc 
show below that even with its limited resolution, numerics 
is a great tool for testing scaling laws. 

Kolmogorov theory provides a scaling law for incom- 
pressible non-magnetized hydrodynamic turbulence (Kol- 
mogorov 1941). This law is true in the statistical sense and 
it provides a relation between the relative velocity vi of fluid 
elements and their separation I, name ly, VI ~ An equiv- 
alent description is to express spectrum E{k) as a function 
of wave number k (~ 1//). The two descriptions are re- 
lated by kE{k) ~ vf . The famous Kolmogorov spectrum 
is E{k) ~ k~^^^. The applications of Kolmogorov theory 
range from engineering research to meteorology (see Monin 
& Yaglom 1975) but its astrophysical applications are poorly 
justified. 

Let us consider incompressible MHD turbulence first. 
There have long been understanding that the MHD tur- 
bulence is anisotropic^ (e.g. Shebalin et al. 1983). A sub- 
stantial progress has been achieved recently by Goldreich & 
Sridhar (1995; hereafter GS95) who made an ingenious pre- 
diction regarding relative motions parallel and perpendicu- 
lar to magnetic field B for incompressible MHD turbulence. 
The GS95 model envisages a Kolmogorov spectrum of ve- 
locity and a scale-dependent anisotropy (see below). These 
relations have been confirmed numerically (Cho & Vishniac 
2000b; Maron &: Goldreich 2001; Cho, Lazarian & Vishniac 
2002a, hereafter CLV02a; see also CLV03a) ; they are in good 
agreement with observed and inferred astrophysical spectra 
(see CLVOSa) . A remarkable fact revealed in CLV02a is that 
fiuid motions perpendicular to B are identical to hydrody- 
namic motions. This provides an essential physical insight 
into why in some respects MHD turbulence and hydrody- 
namic turbulence are similar, while in other respects they 
are diflferent. 

The GS95 model considered incompressible MHD, but 
the real ISM is highly compressible. Literature on the prop- 
erties of compressible MHD is very rich (see CLV03a) . Back 
in 80's Higdon (1984) theoretically studied density fiuctu- 
ations in the interstellar MHD turbulence. Matthacus & 
Brown (1988) studied nearly incompressible MHD at low 
Mach number and Zank & Matthaeus (1993) extended it. 
In an important paper Matthaeus et al. (1996) numerically 



It is not possible to cite all the important papers in the area of 
MHD turbulence. An incomplete list of the references in a recent 
review on the statistics of MHD turbulence by Cho, Lazarian & 
Vishniac (2003a; henceforth CLVOSa) includes about two hundred 
entries. 



explored anisotropy of compressible MHD turbulence. How- 
ever, those papers do not provide urnversal scalings of the 
GS95 type. 

Is it feasible to obtain scaling relations for the compress- 
ible MHD turbulence? Some hints about effects of compress- 
ibility can be inferred from the Goldreich & Sridhar (GS95) 
seminal paper. A more focused discussion was presented in 
the Lithwick & Goldreich (2001) paper which deals with 
electron density fluctuations in the gas pressure dominated 
plasma, i.e. in high 13 regime (/3 = Pgas/Pmag > 1). Incom- 
pressible regime formally corresponds to /3 ^ oo and there- 
fore it is natural to expect that iov (i^ 1 the GS95 picture 
would persist. Lithwick & Goldreich (2001) also speculated 
that for low (3 plasmas the GS95 scaling of slow modes may 
be applicable. An important study of MHD modes in com- 
pressible low /3 plasmas is given in Cho & Lazarian (2002; 
hereafter CL02) where we developed and tested our tech- 
nique of separating different MHD modes. 

In this work, we provide a detailed study of mode 
coupling and scalings of compressible (fast and slow) and 
Alfvenic modes in high j3, intermediate, and low (3 plas- 
mas. Our approach is complementary to that employed in 
direct numerical simulations of astrophysical turbulence. In 
such simulations, e.g. in those dealing with the interstel- 
lar medium (see Vazquez- Semadeni et al. 2000), siirmlations 
of particular astrophysical objects, e.g. molecular clouds, 
are attempted. These simulations provide synthetic maps 
that can be compared with observations. Our goal is to ob- 
tain scaling laws that can also be compared with observa- 
tions. In §2, we describe our approach to the problem in- 
cluding both simple theoretical considerations/expectations 
that motivate our study and the numerical technique that we 
employ. In §3, we present velocity spectra and anisotropics 
for high and low (3 plasmas. In §4, we discuss scalings of 
density and magnetic field. In §5, we present the study of 
viscosity-damped regime of MHD turbulence in compressible 
fluid. This study extends our earlier work (Cho, Lazarian, 
& Vishniac 2002b, henceforth CLV02b) where this regime 
was reported for incompressible flows. In §6, we discuss as- 
trophysical implications of our results, including the rate of 
MHD turbulence decay, relation between the decay rate and 
compressibility, correlation of density and magnetic field, 
and formation of density structures at AU scale. The sum- 
mary is given in §7. 



2 OUR APPROACH 

2.1 Theoretical Considerations 

Let us start with a discussion why isotropic Kolmogorov tur- 
bulence cannot be applicable for describing strongly magne- 
tized gas. Assume that, at some large scale L, the magnetic 
energy and kinetic energy are equal: pVl/2 ~ B^/(47r). Ac- 
cording to the Kolmogorov theory, the kinetic energy at scale 
/ < L is pV^ /2 ~ {l/Lf''-\pVl/2), which is smaller than the 
large scale kinetic energy by a factor of {l/Vf^^. But, mag- 
netic energy density does not diminish as the scale reduces. 
Therefore, at scales smaller than L, hydrodynamic motions 
will not be able to bend magnetic field lines substantially. 

An important observation that leads to understand- 
ing of the GS95 scaling is that magnetic field cannot pre- 
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Figure 1. Separation method. We separate Alfven, slow, and fast modes in Fourier space by projecting the velocity Fourier component 
Vk onto bases ^^i, ^s, and ^f, respectively. Note that g^i = —<fi. Slow basis ^3 and fast basis lie in the plane defined by Bq and k. 
Slow basis lies between —8 and k||. Fast basis lies between k and kj^. 



vent mixing motions of magnetic field lines if the motions 
are perpendicular to the magnetic field. Those motions will 
cause, however, waves that will propagate along magnetic 
field lines. If that is the case, the time scale of the wave-like 
motions, i.e. ~ /||/Va, where is the characteristic size of 
the perturbation along the magnetic field and Va = B/y^iirp 
is the local Alfven speed, will be equal to the hydrody- 
namic time-scale, l±/vi, where l± is the characteristic size 
of the perturbation perpendicular to the magnetic field. 

The mixing motions are hydrodynamic-like'^ and therefore 

1/3 

obey Kolmogorov scaling vi cc . Combining the two re- 
lations above we can get the GS95 anisotropy, l\\ oc f"/^^ (or 
fc|| oc k^/^^ in terms of wave-numbers). If we interpret l\\ as 
the eddy size in the direction of the local* magnetic field and 
l± as that in the perpendicular directions, the relation im- 
plies that smaller eddies are more elongated (see Appendix 
B for illustration of scale-dependent anisotropy). 

How is this idealized incompressible model related to 
the actual, e.g. interstellar, turbulence? Compressible MHD 
turbulence is a highly non-linear phenomenon and it has 
been thought that different types of perturbations or modes 
(Alfven, slow and fast) in compressible media are strongly 
coupled. Nevertheless, one may question whether this is true. 
A remarkable feature of the GS95 model is that Alfven per- 
turbations cascade to small scales over just one wave period, 
while the other non-linear interactions require more time. 
Therefore one might expect that the non-linear interactions 
with other types of waves should affect Alfvenic cascade only 
marginally. Moreover, as the Alfven waves are incompress- 



^ Recent simulations (Cho et al. 2003) suggest that perpendicu- 
lar mixing is indeed efficient for mean magnetic fields of up to the 
equipartition value. This corresponds to our earlier result that 
high order velocity statistics of MHD turbulence in the perpen- 
dicular directions is very similar to that of hydrodynamic one 
(CLV02a). 

* The concept of local is crucial. The GS95 scalings are obtained 
only in the local frame of magnetic field, as this is the frame where 
magnetic field are allowed to be mixed without being opposed by 
magnetic tension. 



ible, the properties of the corresponding cascade may not 
depend on the sonic Mach number. 

The generation of compressible motions (i.e. rarfiai com- 
ponents in Fourier space) from Alfvenic turbulence is a 
measure of mode coupling. How much energy in compress- 
ible motions is drained from Alfvenic cascade? According 
to the closure calculations (Bertoglio, Bataille, & Marion 
2001; see also Zank & Matthaeus 1993), the energy in com- 
pressible modes in hydrodynamic turbulence scales as ~ 
if Ms < 1. We may conjecture that this relation can be 
extended to MHD turbulence if, instead of Mf, we use 
~ {5V)\/{a^ + Vl). (Hereinafter, we define Va = Bq/^/^, 
where Bo is the mean magnetic field strength.) However, as 
the Alfven modes are anisotropic, this formula may require 
an additional factor. The compressible modes are generated 
inside the so-called Goldreich-Sridhar cone, which takes up 
~ {5V)a/Va of the wave vector space. The ratio of com- 
pressible to Alfvenic energy inside this cone is the ratio given 
above. If the generated fast modes become isotropic (see 
below), the diffusion or, "isotropization" of the fast wave 
energy in the wave vector space increase their energy by a 
factor of ~ Va/{5V)a- This results in 



Vl + {&V)a 
Va 



(1) 



where {5V)l^d and {5V)\ are energy of compressible^ and 
Alfven modes, respectively. Eq. suggests that the drain 
of energy from Alfvenic cascade is marginal when the am- 
plitudes of perturbations are weak, i.e. {5V)a <Si Va- 

If Alfven cascade evolves on its own, it is natural to 
assume that slow modes exhibit the GS95 scaling. Indeed, 
slow modes in gas pressure dominated environment (high /3 
plasmas) are similar to the pseudo- Alfven modes in incom- 
pressible regime (see GS95; Lithwick & Goldreich 2001). The 
latter modes do follow the GS95 scaling. In magnetic pres- 
sure dominated environments (low 13 plasmas), slow modes 



It is possible to show that the compressible modes inside the 
Goldreich-Sridhar cone are basically fast modes. 
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are density perturbations propagating with the sound speed 
a parallel to the mean magnetic field (see equation (IA31II '). 
Those perturbations are essentially static for a ^ Va- 
Therefore Alfvenic turbulence is expected to mix density 
perturbations as if they were passive scalar. This also in- 
duces GS95 spectrum. 

The fast waves in low /3 regime propagate at Va irre- 
spectively of the magnetic field direction. In high /3 regime, 
the properties of fast modes are similar, but the propaga- 
tion speed is the sound speed a. Thus the mixing motions 
induced by Alfven waves should marginally affect the fast 
wave cascade. It is expected to be analogous to the acoustic 
wave cascade and hence be isotropic. 

For most part of this paper, we shall assumed that 
SV ~ SB/y^inp ~ Bo /v'47rp, where SB is the r.m.s. strength 
of the random magnetic field. This is less restrictive than it 
might appear, since as long as there is some scale I' in the 
turbulent cascade where u;/ ~ B/y^inp we can take L = I' , 
Vl = vii and use this model of turbulence for all smaller 
scales. Suppose that we have a turbulent system initially 
threaded by a mean magnetic field only. If initially the tur- 
bulent energy is larger than the magnetic energy of the mean 
field, we are in the regime of so-called superAlfvenic tur- 
bulence. In this regime the growth of the magnetic field is 
expected through so called "turbulent dynamo" (see Batch- 
elor 1950; Brandenburg et al. 1996; Cho & Vishniac 2000a). 
Initially, the growth of magnetic energy is most active at the 
scale roughly an order of magnitude larger than the dissi- 
pation scale and the magnetic spectrum peaks at the scale. 
As the magnetic energy grows, the magnetic back reaction 
becomes important at the scale and the peak of the mag- 
netic power spectrum moves to larger scales. Finally, when 
equipartition between the kinetic and magnetic energy den- 
sities occurs at a scale somewhat (factor of 2 or 3) smaller 
than the kinetic energy peak, turbulence reaches a statis- 
tically stationary state. This agrees well with the results 
of incompressible simulations in Cho & Vishniac (2000a). 
However, this is not universally accepted idea. For example, 
Padoan & Nordlund (1999) reported that, in their compress- 
ible simulations, the power spectrum of p^^'^v has a shal- 
lower slope then those of velocity and magnetic field, which 
implies that equipartition between kinetic magnetic energy 
cannot be reached at small scales. At scales smaller than the 
equipartition scale, the turbulence becomes subAlfvenic and 
our earlier considerations should be applicable. 

The arguments above suggest that compressible MHD 
turbulence should demonstrate well defined scaling rela- 
tions. Below we test those arguments and reveal scaling re- 
lations for different /3s and Mach numbers. 



2.2 Numerical scheme 

We use a third-order accurate hybrid essentially non- 
oscillatory (ENO) scheme (see CL02) to solve the ideal 
isothermal MHD equations in a periodic box: 

dp/dt + V ■ (pv) = 0, (2) 
d^/dt V ■ Vv + p~^V{ap) - (V X B) X B/47rp = f, (3) 

9B/at- V X (v X B) = 0, (4) 

with V • B = and an isothermal equation of state. Here f 
is a random large-scale driving force, p is density, v is the 



velocity, and B is magnetic field. The rms velocity 5V is 
maintained to be approximately unity (in fact 5V ~ 0.7), so 
that V can be viewed as the velocity measured in units of 
the r.m.s. velocity of the system and B/-^47rp as the Alfven 
velocity in the same units. The time t is in units of the 
large eddy turnover time (~ L/5V) and the length in units 
of L, the scale of the energy injection. The magnetic field 
consists of the uniform background field and a fiuctuating 
field: B = Bo -f b. 

For mode coupling studies (Fig. |5J , we do not drive tur- 
bulence. For scaling studies, we drive turbulence solenoidally 
in Fourier space and use 216^ points, Va = Bo /\/47rp = 1, 
and po = 1. The average rms velocity in statistically sta- 
tionary state is 5V ~ 0.7. 

For our calculations we assume that Bo/y^inp ~ 
SB/^4:Tvp ~ SV. In this case, the sound speed is the control- 
ling parameter and basically two regimes can exist: super- 
sonic and subsonic. Note that supersonic means low-beta 
and subsonic means high-beta. When supersonic, we con- 
sider mildly supersonic (or, mildly low-/3) and highly super- 
sonic (or, very low-/3)®. 

2.3 Separation of MHD modes 

Three types of waves exist (Alfven, slow and fast) in com- 
pressible magnetized plasma. The slow, fast, and Alfven 
bases that denote the direction of displacement vectors for 
each mode are given by 

(5) 
(6) 
(7) 



oc (-l-Fa-VD)fc||k| 
if oc {-1 + a + \/D)k\\k\ 



-f (1 + Q - VD)fcxkx 
+ {l + a + VD)k_Lk±^ 
^A = -ifi = 'k± X kii 



where D = {1 + a)^ 4acos9, a = /V^ = /3(7/2), 9 is 
the angle between k and Bo, and ip is the azimuthal basis in 
the spherical polar coordinate system. These are equivalent 
to the expression in CL02: 



oc fc||k|| -f 



D~f3/2 



I k I 



D + /3/2 



fciikii -I- fcxki 



(8) 



(9) 



(Note that 7=1 for isothermal case.) 

Slow and fast velocity components can be obtained by 
projecting velocity Fourier component Vk onto and £,f, 
respectively. In Appendix A, we discuss how to separate slow 
and fast magnetic modes. We obtain energy spectra using 
this projection method. When we calculate the structure 
functions (e.g. for Alfven velocity) in real space, we first 
obtain the Fourier components using the projection and, 
then, we obtain the real space values by performing Fourier 
transform. 



° The terms "mildly" and "highly" are somewhat arbitrary 
terms. We consider these two supersonic cases to cover a broad 
range of parameter space. Note that Boldyrev, Nordlund, & 
Padoan (2002b) recently provided a Mach number dependence 
study of the compressible MHD turbulence statistics where only 
two regimes are manifest: essentially incompressible and essen- 
tially compressible shock-dominated (with smooth transition at 
some Ms of order unity). 
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Figure 2. Mode coupling studies. {a)left: Square of the r.m.s. velocity of the compressible modes. We use 144^ grid points. Only Alfvcn 
modes are allowed as the initial condition. "Pluses" are for low cases (0.02 ^ /3 ^ 0.4). "Diamonds" are for high /3 cases (1^/3^ 20). 
{b)middle: Generation of fast modes. Snapshot is taken at t=0.06 from a simulation (with 144^ grid points) that started off with Alfven 
modes only. Initially, j3 (ratio of gas to magnetic pressure, Pg/ Pmag) = 0.2 and Mg (sonic Mach number) ~ 1.6. {c)right: Comparison 
of decay rates. Decay of Alfven modes is not much affected by other (slow and fast) modes. We use 216^ grid points. Initially, /3 = 0.02 
and Ms ~ 4.5 for the solid line and Ms ~ 7 for the dotted line. Note that initial data axe, in some sense, identical for the solid and the 
dotted lines. The sonic Mach number for the solid line is smaller because we removed fast and slow modes from the initial data before 
the decay simulation. For the dotted line, we did not remove any modes from the initial data. 
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Figure 3. Low /3 (/3 r^ 0.2 and Ms ~2.3). Scalings relations. Results from driven turbulence with Ma (Alfven Mach number) r~ 0.7, 
and 216^ grid points. (Va = Bo / V^Trp = 1; a (sound speed) = VO.l; &V ^ 0.7.) (a.) Upper-left: Spectra of Alfven modes follow a 
Kolmogorov-like power law. (b) Middle-left: The second-order structure function {SF2) for velocity of jVlfvcu modes shows anisotropy 
2/3 2/3 

similar to the GS95 (r|| oc Vj^ or fcy oc kj^ ). The structure functions are measured in directions perpendicular or parallel to the local 
mean magnetic field in real space. We obtain real-space velocity and magnetic fields by inverse Fourier transform of the projected fields. 
(c) Lower-left: SF2 on the parallel axis and on perpendicular axis for Alfven modes velocity, (d) Upper-middle: Spectra of slow modes also 
follow a Kolmogorov-like power law. {e) Middle-middle: Slow mode velocity shows anisotropy similar to the GS95. (i) Lower-middle: SF2 

on the parallel axis and on perpendicular axis for slow modes velocity, (g) Upper-right: Spectra of fast modes are compatible with the IK 
spectrum, (h) Upper-middle: The SF2 of fast modes velocity shows isotropy. Fast mode magnetic field also shows isotropy. (i) Lower-right: 
SF2 on the parallel axis and on perpendicular axis for fast modes velocity. 
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Figure 4. High /3 (/3 --^ 4 and Ma --^0.35). Va = Bq/'JAtxp = I. a (sound speed) = \/2. 5V ~ 0.7. See caption in Fig. 13 Alfven and slow 
modes follow the GS95 scalings. Fast modes are isotropic. 



3 VELOCITY SCALINGS 
3.1 Mode coupling 

In CL02 we demonstrated the decoupling of Alfven and fast 
modes in low /3 plasmas. Here we substantially extend the 
CL02 analysis. As mentioned above, the coupling of com- 
pressible and incompressible modes is crucial. If Alfvenic 
modes produce a copious amount of compressible modes, 
the whole picture of independent Alfvenic turbulence fails. 
However, our calculations show that the amount of energy 
drained into compressible motions is negligible, provided 
that either the external magnetic field or the gas pressure 
is sufficiently high. Fig. |2K suggests that the generation of 
compressible motions follows equation Fast modes also 
follow a similar scaling, although the scatter is a bit larger. 
The marginal generation of compressible modes is in agree- 
ment with earlier studies by Boldyrev et al. (2002b) and 
Porter, Pouquet, & Woodward (2002), where the velocity 
was decomposed into a potential component and a solenoidal 
component. See Fig. |2ji for the values of ~ x (=the ra- 
tio of the mean square potential to solenoidal velocity). 
Fig. I^t demonstrates that fast modes are initially generated 
anisotropically, which supports our theoretical consideration 
in 

Fig. |5J; shows that dynamics of Alfven modes is not 
affected by slow modes. We first perform a driven turbulence 



simulation with 216"^ grid point, (3 ~ 0.02, and Ms ~ 7. 
Then, after it has reached a statistically stationary state, we 
stop the run and save the data. Using these data, we perform 
two decay simulations. For one (the solid line), we remove 
all slow and fast modes and let the turbulence decay. For the 
other (the dotted line), without removing any modes, we just 
let the turbulence decay. The solid line in the figure is the 
energy in Alfven modes when we start the decay simulation 
with Alfven modes only. The dotted line is the Alfven energy 
when we start the simulation with all modes. This result 
confirms that Alfven modes cascade is almost independent of 
slow and fast modes. In this sense, coupling between Alfven 
and other modes is weak. 



3.2 High-/? and mildly supersonic low-/3 regimes 

Alfven Modes. — Fig. (low-/3) and Fig. ^ (high-/3) 
show that the power spectra of Alfven waves follow a Kol- 
mogorov spectrum: 

Spectrum of Alfven Waves: E^{k) (x kj_^^^ . (10) 

In Fig. {middle-left panel) and Fig.jljj {middle-left panel), 
we plot contours of equal second-order structure function 
for velocity {SF2{r) =< |v(x-|-r) — v(x)|^ >avg. over x) ob- 
tained in local coordinate systems in which the parallel axis 
is aligned with the local mean field (see Cho & Vishniac 
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2000b; CLV02a; Maron & Goldreich 2001). The SF2 along 
the axis perpendicular to the local mean magnetic field fol- 
lows a scahng compatible with r'^^^ . The SF2 along the axis 
parallel to the local mean field follows steeper scaling 
(Fig. I2J; and Fig. lower-left panels). The results show 
reasonable agreement with the GS95 model for incompress- 
ible MHD turbulence, 

Amsotropy of Alfven Waves: ry oc r^/^ , or fcy oc , (11) 

where ry and r± are the semi-major axis and semi-minor 
axis of eddies, respectively (Cho & Vishniac 2000b). 

Slow waves. — The incompressible limit of slow waves is 
pseudo- Alfven waves. Goldreich & Sridhar (1997) argued 
that the pseudo- Alfven waves are slaved to the shear- Alfven 
(i.e. ordinary Alfven) waves, which means that pseudo- 
Alfven modes do not cascade energy for themselves. Lith- 
wick & Goldreich (2001) made similar theoretical arguments 
for high 13 plasmas and conjectured similar behaviors of slow 
modes in low /3 plasmas. We confirmed that similar argu- 
ments are also applicable to slow waves in low l3 plasmas 
(CL02). Indeed, power spectra in Fig.j^i and Fig.|l]l {upper- 
middle panels) are consistent with: 

Spectrum of Slow Modes: E'' (k) (x k~^^^ . (12) 

In Fig. 1^ and Fig. 2t {middle-middle panels), contours of 
equal second-order velocity structure function {SF2), rep- 
resenting eddy shapes, show scale-dependent anisotropy: 
smaller eddies are more elongated. The results show rea- 
sonable agreement with the GS95-type anisotropy 

Amsotropy of Slow Modes: fcy oc k^/'^ , or ry oc r^'', (13) 

where ry and r± are the semi-major axis and semi-minor 
axis of eddies, respectively. 

Fast waves. — Fig. I^i and Fig. 2h {middle-right panels) 
show fast modes are isotropic. The resonance conditions for 
the interacting fast waves are lu\-\-lu2 ~ i^a and ki -I- k2 = 
ks. Since uo oc k for the fast modes, the resonance condi- 
tions can be met only when all three k vectors are coUinear. 
This means that the direction of energy cascade is radial in 
Fourier space. This is very similar to acoustic turbulence, 
turbulence caused by interacting sound waves (Zakharov 
1967; Zakharov & Sagdeev 1970; L'vov, L'vov, & Pomyalov 
2000). Zakharov & Sagdeev (1970) found E{k) oc fc'^/^ 
However, there is debate about the exact scaling of acous- 
tic turbulence. Here we cautiously claim that our numerical 
results are compatible with the Zakharov & Sagdeev scaling: 

Spectrum of Fast Modes: {k) ^ k'''^'^ . (14) 

The eddies are isotropic (see also Appendix B). 

3.3 Highly supersonic low-/3 case 

The results for low-/3 in the previous subsection are for a 
Mach number of ~ 2.3. In this subsection, we present results 
for a Mach number of ~ 7. Obviously shock formation is 
faster when the Mach number of the system is high. We also 
expect that turbulent motions can compress/disperse the 
gas more easily when the Mach number is high. As a result, 
we expect higher density fiuctuations when Mach number is 



higher. Thus we check the scaling relations for high Mach 
number fiuids. 

Fig. Oshows that most of the scaling relations that hold 
true in mildly supersonic flows are still valid in the highly 
supersonic case. Especially anisotropy of Alfven, slow, and 
fast modes is almost identical to the one in the previous 
section. However, the power spectra for slow modes do not 
show the Kolmogorov slope. The slope is close to —2, which 
is suggestive of shock formation. At this moment, it is not 
clear whether or not the —2 slope is the true slope. In other 
words, the observed —2 slope might be due to the limited 
numerical resolution. Runs with higher numerical resolution 
should give the definite answer. 



3.4 SuperAlfvenic turbulence 

So far we considered "sub-Alfvenic" turbulence in which 
the Alfven speed associated with the mean magnetic field 
is slightly faster than the r.m.s. fiuid velocity. If the op- 
posite case is true (i.e. if the mean field Bq is weak), 
the turbulence is called "super- Alfvenic" . In super- Alfvenic 
turbulence, large scale magnetic field lines can show very 
chaotic structures. Whether or not the ISM turbulence is 
sub-Alfvenic is still a controversial issue. 

We mentioned in H2. II that, even in the case of super- 
Alfvenic turbulence, we can find some scale I' in the tur- 
bulent cascade where vii ~ B/y/Aup and we can apply our 
model of sub-Alfvenic turbulence for all smaller scales. Fig.|S| 
supports this idea. The contours in the figure are the sec- 
ond order structure functions (SF2) of velocity and magnetic 
field. We do not use mode decomposition. Nevertheless, the 
velocity SF2 refiects scalings of Alfven and slow modes, be- 
cause fast modes are weaker than Alfven and slow modes. 
As expected, anisotropy emerges at small scales. This is very 
similar to incompressible case (Cho & Vishniac 2000a). 

Fig- EI a) shows power spectra. We notice that the power 
spectra of velocity and magnetic field have different shapes. 
The velocity power spectrum is larger than the magnetic 
one near the energy injection scale at fc ~ 2.5. However, for 
larger k's {k > 6), the magnetic power spectrum is larger 
than the velocity one. This behavior is well known in in- 
compressible simulations with unit magnetic Prandtl num- 
ber (see, for example, Kida, Yanase, & Mizushima 1991; 
Cho & Vishniac 2000a). A similar behavior is also observed 
in earlier compressible simulations (see, for example, Bran- 
denburg et al. 1996). Cho & Vishniac (2000a) argued that 
the transition from Ev{k) > Eb{k) to Ev{k) < Eb{k) occurs 
at a wavelength 2 or 3 time larger than that of the energy 
injection scale. 

A careful look at Fig. ^a,) reveals that, for fc < 8, the 
power spectrum of velocity declines faster than Kolmogorov 
as k increases. This is not very surprising because magnetic 
field has more power than velocity at large k's and, therefore, 
can affect the velocity power spectrum at small k's. Kida et 
al. (1991) claimed that the sum of Ev{k) + Ei,{k) roughly 
follows Kolmogorov spectrum in their incompressible sim- 
ulations. If this is true, then the velocity power spectrum 
should have a spectrum steeper than Kolmogorov at small 
k's because many simulations have shown that the magnetic 
power spectrum is significantly fiatter than the velocity one 
at small k's when the mean field Bo is weak (Kida et al. 
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Figure 5. Highly supersonic low /3 (/3 ~ 0.02 and Ms ^7). Va = Bo / inp = I. a (sound speed) = 0.1. 5V ^ 0.7. Alfven modes follow 
the GS95 scalings. Slow modes follow the GS95 anisotropy. But velocity spectrum of slow modes is uncertain. Fast modes arc isotropic. 




1991; Brandenburg et al. 1996; Cho & Vishniac 2000a). Af- larger scales than ours and their sonic Mach number is larger 
ter ~ 8, it seems that the velocity power spectrum gets than ours. Further parameter study is absolutely necessary, 
flatter. 



Boldyrev, Nordlund, & Padoan (2002a) also obtained 
velocity power spectrum steeper than Kolmogorov in 
their supersonic super-Alfvenic MHD simulations. They at- 
tributed the steep spectrum to different intermittency prop- 
erties compared to the incompressible case. Since they used 
a different simulation set-up, we do not directly compare our 
results and theirs. For example, their turbulence is driven at 



3.5 How good is our technique? 

The technique described in H2..'-{l is statistical in nature. That 
is, we separate each MHD mode with respect to the mean 
magnetic field Bq. This procedure is affected by the wan- 
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Figure 7. Comparison between Fourier space metliod and real space method. {a,)left: From real space calculation. Ms ~ 7. (h)middle: 
Solid: Fourier space. Dashed: real space. Ms ~ 7. {c)right: Similar plot for Ms ~ 2.3. 



dering of large scale magnetic field lines, as well as density 
inhomogeneities^ . 

Nevertheless, we can show that our technique gives sta- 
tistically correct results. In low /3 regime, the velocity of a 
slow mode is nearly parallel to the local mean magnetic field. 
Therefore, for low /3 plasmas, we can obtain velocity statis- 
tics for slow modes in real space as follows. First, we iden- 
tify the direction of the local mean magnetic field using the 
method described in Cho & Vishniac (2000b). Second, we 
calculate the second order structure function for slow modes 
by the formula vSF2(r) =< ] (v(x + r) - v(x)) ■ B;]^ >, 
where B; is the unit vector along the local mean field. 

Fig. |7[a) shows the contours obtained by the method 
for the high Mach number run. In Fig.|7[b), we compare the 
result obtained this way (dashed lines) and our technique 
described in H2.3I (solid lines). We also show a similar plot 
for the mildly supersonic case in Fig. |7[c). These results 
confirm that the method described in H2.3l gives statistically 
correct scaling relations. 



4 LINEAR ESTIMATES OF DENSITY AND 
MAGNETIC FIELD FLUCTUATIONS 

In this section, we estimate the magnitude of the r.m.s. 
fluctuations of magnetic field and density. Figures El El 
and 1^ show that magnetic field and density have spectral 
indexes similar to those of velocity. We also expect that 
isotropy/anisotropy of magnetic field is similar to that of 
velocity. Therefore, we do not discuss these quantities here. 
However, anisotropy of density shows different behavior. Fig. 
|51 shows structure of density. Density shows anisotropy for 
the high /3 case. But, for low (3 cases, density shows more or 
less isotropic structures. We suspect that shock formation is 
responsible for the isotropization of density. 

To estimate the r.m.s. fluctuations, we use the following 
linearized continuity and induction equations: 



\pk\ = {poVk/c)\k-(,\, 

\hk\ = (Bot'fc/c)lB„ X CI, 



(15) 
(16) 



^ One way to remove the effect by the wandering of field lines is 
to drive turbulence anisotropically in such a way as fcj^ i^&V ~ 
A;|| ^Va, where kj^ j^ and stand for the wavelengths of the 

driving scale and 5V is the r.m.s. velocity. By increasing the 
l'^,L/k\\,L ratio, we can reduce the degree of mixing of different 
wave modes. 



where c denotes propagation speed of slow or fast wave 
(equation \A22V \. From this, we obtain the r.m.s. fluctu- 
ations 



{6p/po)f 
{5B/Bo)s 
{5B/Bo)f 



{SV)s{\k-^s/Cs\), 

{5V)f{\k-if/cf\}, 
{SV)s{\Bo X is/cs\), 
{5V)f{\Bo X if/cf\), 



(17) 
(18) 
(19) 
(20) 



where angled brackets denote a proper Fourier space aver- 
age. Generation of slow and fast modes velocity {{SV)s and 
{SV)f) depends on driving force. Therefore, we may simply 
assume that 



{5V)a ~ {SV)s ~ {SV)f, 



(21) 



where we ignore constants of order unity. However, when 
we consider mostly incompressible driving, the generation 
of fast modes follows equation Q. In this case, the am- 
plitude of fast mode velocity is reduced by a factor of 



(SV) 



Va 



-1/2 



(5V)a ~ {5V)s 



Vl + a"" {5V)a 



Va 



1/2 



(22) 



When we assume {SV)a ~ Va, equation 1221 reduces to 
equation H21|l in low /3 plasmas. 



Va- Using equations 



4.1 Low-/3 case 

In this limit, Cs ~ acos9 and c/ 
HASH and 1A32II . we obtain 

{5p/po)s ~ {5V) s{\ cos 6 / cs\) {5V)s/ a, (23) 

{Sp/po)f = {SV)f{\sme/cf\)^{SV)f/VA, (24) 

{SB/Bo)s = (<51/),(|acosesin6'/c^j) ~ a((5V)^/a, (25) 

{5B/Bo)f = {SV)f{\l/cf\)^{5V)f/VA, (26) 

where we ignore cos6''s or sink's. 

When we assume {5V)a ~ {SV)s 



{SV)j 



{Sp/po)s 

{Sp/po)f 
{5B/Bo)s 
{SB/Bo)f 



Ms, 

y/pMs, 

l3Ms, 

^/pMs 



Va, we get 
(27) 
(28) 
(29) 
(30) 



Therefore, in low 13 plasmas, slow modes give rise to most of 
density fluctuations (CLQ2). On the other hand, magnetic 
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Figure 8. Density structures. {a,)left: Ms ~ 0.35 (high f3). {h)middle: Ms ~ 2.3 (mildly supersonic low /3). (c)right: Ms ~ 7 (highly 
supersonic low (3). 



fluctuation by slow modes is smaller than that by fast modes 
by a factor of ^/p. 



4.2 High-/? case 

In this limit, Cs ~ Vacos6' and Cf 
<A33l l and 1A34II . we obtain 



a. Using equations 



[&p/po)s ~ (5V)»{|cos6lsin6l/(ac,)|) 

~ {VA/a)[SV)s/a, 

{5p/po)f = {SV)f{\l/cf\)r^{SV)f/a, 

{5B/Bo)s = {SV)s{\cose/cs\) {5V)s/Va, 

{5B/Bo)f = {SV)f{\sme/cf\)^{SV)f/a, 

where we ignore cos6''s or sink's. 

Let us just assume that {5V)a ~ {SV)s 
Mr^{5V)f (cf. equation Then we have 



(31) 
(32) 
(33) 
(34) 



Va 



{Sp/po)f 
{SB/Bo)s 
{5B/Bo)f 



Ms/^/^' 
Ml 
Oil), 
Ml 



Mi 



(35) 
(36) 
(37) 
(38) 



The density fluctuation associated with slow modes is ~ Ml 
when {5V)a ~ {5V)a ~ Va- This is consistent with Zank 
& Matthaeus (1993). The ratio of {5p)s to {5p)f is of or- 
der unity. Therefore, both slow and fast modes give rise to 
similar amount of density fluctuations. Note that this ar- 
gument is of order-of-magnitude in nature. In fact, in our 
simulations for the high (5 case, the r.m.s. density fluctua- 
tion by slow modes is about twice as large as that by fast 
modes. When we use equation I21II . we have a different re- 
sult: {5p)s ~ {Va/ a){5p) f < {5p)f. It is obvious that slow 
modes dominate magnetic fluctuations: {SB)s > {SB)f for 
both equations 1211 1 and 12211 . 



the end of the hydrodynamic cascade, but in MHD turbu- 
lence this is not the end of magnetic structure evolution. 
For viscosity much larger than resistivity, v ^ rj, there 
will be a broad range of scales where viscosity is impor- 
tant but resistivity is not. On these scales magnetic field 
structures will be created by the shear from non-damped 
turbulent motions, which amounts essentially to the shear 
from the smallest undamped scales. The created magnetic 
structures would evolve through generating small scale mo- 
tions. As a result, we expect a power-law tail in the mag- 
netic energy distribution, rather than an exponential cut- 
off. Cho, Lazarian, & Vishniac (CLV02b) performed numer- 
ical simulations of turbulence in this regime threaded by a 
strong {Bo/^/4np ~ SV) mean magnetic field and reported 
that this regime possesses completely different scaling re- 
lations and anisotropic structures compared with ordinary 
MHD turbulence. Further research showed that there is a 
smooth connection between this regime and small scale tur- 
bulent dynamo in high magnetic Prandtl number fluids (see 
Schekochihin et al. 2002)*. 

In partially ionized gas neutrals produce viscous damp- 
ing of turbulent motions. In the Cold Neutral Medium (see 
Draine & Lazarian 1999 for a list of the idealized phases) 
this produces damping on the scale of a fraction of a par- 
sec. The magnetic diffusion in those circumstances is still 
negligible and exerts an influence only at the much smaller 
scales, ~ lOOfcm. Therefore, there is a large range of scales 
where the physics of the turbulent cascade is very different 
from the conventional MHD turbulence picture. 

CLV02b explored this regime numerically with a grid 
of 384"^ and a physical viscosity for velocity damping. The 
kinetic Reynolds number was around 100. We achieved a 
very small magnetic diffusivity by the use of hyper-diffusion. 
The result is presented in Fig. A theoretical model for 
this regime and its consequences for stochastic reconnection 
(Lazarian & Vishniac 1999) can be found in Lazarian, Vish- 
niac, & Cho (2003). It explains the spectrum E{k) ~ as 



5 VISCOSITY-DAMPED REGIME OF MHD 
TURBULENCE 

In hydrodynamic turbulence viscosity sets a minimal scale 
for motion, with an exponential suppression of motion on 
smaller scales. Below the viscous cutoff the kinetic energy 
contained in a wavenumber band is dissipated at that scale, 
instead of being transferred to smaller scales. This means 



" It is worth noting that, motivated by the analogy between time 
evolution equations for magnetic field and for vorticity, Batchelor 
(1950) first argued that small magnetic fields can be amplified 
when viscosity is larger than magnetic diffusion (i.e. magnetic 
Prandtl number > 1). Although the analogy was later proved to 
be physically wrong, the high magnetic Prandtl number dynamo 
has been studied by many researchers (e.g. Kulsrud & Anderson 
1992; Kinney et al. 2000). 
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Figure 9. Viscous damped regime (viscosity > magnetic difTusivity ) . Due to large viscosity, velocity damps after k ~ 10. (a) Left: 
Incompressible case with 384'' grid points. Magnetic spectra show a shallower slope {Ei,{k) oc k~^) below the velocity damping scale. 
We achieve a very small magnetic difTusivity through the use hyper-diffusion. From Cho, Lazarian, & Vishniac (2002b). (b) Right: 
Compressible case with 216'^ grid points. Magnetic and density spectra show structures below the velocity damping scale at fc ~ 10. The 
structures are less obvious than the incompressible case because it is relatively hard to achieve very small magnetic diffusivity in the 
compressible run. 



a cascade of magnetic energy to small scales under the influ- 
ence of shear at the marginally damped scales. The spectrum 
is similar to that of the viscous-convective range of passive 
scalar in hydrodynamic turbulence (see, for example, Batch- 
elor 1959; Lesieur 1990), although the study in Lazarian, 
Vishniac, & Cho (2003) suggests that the physical origin of 
it is different. A study confirming that the spectrum 
is not a bottleneck effect is presented in Cho, Lazarian, & 
Vishniac (2003b). The mechanism is based on the solenoidal 
motions and therefore the compressibility should not alter 
the physics of this regime of turbulence. 

We show our results for the compressible fluid in Fig|^. 
We use the same physical viscosity as in incompressible case 
(see CLV02b). We rely on numerical diffusion, which is much 
smaller than physical viscosity, for magnetic field. The iner- 
tial range is much smaller due to numerical reasons, but it is 
clear that the viscosity-damped regime of MHD turbulence 
persists. The magnetic fluctuations, however, compress the 
gas and thus cause fluctuations in density. This is a new 
(although expected) phenomenon compared to our earlier 
incompressible calculations. These density fluctuations may 
have important consequences for the small scale structure of 
the ISM. 



6 ASTROPHYSICAL IMPLICATIONS OF OUR 
RESULTS 

6.1 Parameter range explored 

Parameters of astrophysical plasmas differ substantially for 
different astrophysical systems: from extremely high (3 to 
extremely low /3. Turbulence in some systems is expected to 
be superAlfvenic, in others it is expected to be subAlfvenic. 
Moreover, there is an ongoing controversy on what to ex- 
pect and where. For instance, while high l3 was considered a 
default for many phases of Milky Way ISM, recent observa- 
tions by Beck (2002) suggest that the plasmas there may be 



low p. Therefore it is essential to have clear understanding 
of MHD turbulence for as large parameter space as possible. 

SuperAlfvenic regime. — We have argued above that 
the difference between superAlfvenic and subAlfvenic turbu- 
lence is not as substantial as it is frequently thought. The 
difference between the two regimes stems from ratio of mag- 
netic field to kinetic energies. However, as we mentioned 
in §2, if kinetic energy density exceed magnetic field energy 
density, the hydromagnetic motions drive turbulent dynamo. 
This enhances magnetic field energy density up to approxi- 
mately equipartition value. Thus the difference between the 
superAlfvenic and subAlfvenic regimes amounts not to the 
energy of the magnetic field, but to the global level of field 
organization, e.g. to the magnetic field reversals etc. Our re- 
sults in H3. 41 suggest that the basic properties of the MHD 
turbulence in the subAlfvenic and superAlfvenic regimes are 
similar. This, however, does not preclude the intermittency 
of MHD turbulence being very different. The latter prop- 
erty can be tested using scalings of the higher order veloc- 
ity correlations SFp(r) = {|v(x + r) - v(x)|'') oc r'^^^K The 
corresponding scaling suggested by She-Leveque (1994) con- 
tains three parameters (see Politano & Pouquet 1995; Miiller 
& Biskamp 2000): g is related to the scaling Svi ~ l^^^, x 
related to the energy cascade rate t^^ ~ , and C, the 
co-dimension of the dissipative structures: 

C(P) = P/gil - :c) + C (1 - (1 - x/Cr/') . (39) 

Miiller & Biskamp (2000) proposed that 3D incompressible 
MHD turbulence for zero mean field has Kolmogorov g and 
X, while the the dissipation happens in the sheet-like struc- 
tures, i.e. C = 3 — 2 = 1. Using eq. I39II they obtained an 
excellent fit for their incompressible data. Later Boldyrev 
(2002) made the same assumption^ that C — 1 (and same g 



^ The physical motivation of Miiller & Biskamp (2000) and 
Boldyrev (2002) for choosing the same value of C are differ- 
ent, however. Boldyrev (2002) identifies the dissipation structures 
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and x) for the compressible turbulence and Boldyrev, Nord- 
lund, & Padoan (2002a) obtained an excellent fit to their 
compressible data. It appears surprising that incompressible 
MHD (Miiller & Biskamp 2000) and supersonic compress- 
ible MHD (Boldyrev et al. 2002a) have similar intermittency 
structures. This issue is discussed in Cho, Lazarian, & Vish- 
niac (2003b). 

We would like to note, however, that our considerations 
about superAlfvenic turbulence are not applicable to the 
transient regimes when superAlfvenic motions are in the 
process of generating magnetic field and magnetic energy 
does not have time to come to a rough equipartition to the 
kinetic energy. Some parts of the ISM can well be in the 
transient regime for which temporary pV^/2 ^ B^/Stt. 

/3 = I case. — We provided theoretical considerations 
for both /9 2> 1 and /3 ^ I cases. What about the inter- 
mediate cases of /3 ~ 1? Would the scaling of modes and 
mode coupling be different? To test this case we performed 
calculations for (3 = 1. The results in Fig. 1101 show that the 
scaling relations that hold true for mildly low /? regime are 
also applicable for /3 = 1 regime. 

6.2 Fundamental questions 

How fast does MHD turbulence decay? This question 
has fundamental implications for star formation (see McKee 
1999). Indeed, it was thought originally that magnetic fields 
would prevent turbulence from fast decay. Later (see Mac 
Low et al. 1998; Stone et al. 1998; and review by Vazquez- 
Semadeni et al. 2000) this was reported not to be true. How- 
ever, fast decay was erroneously associated with the coupling 
between compressible and incompressible modes. The idea 
was that incompressible motions quickly transfer their en- 
ergy to the compressible modes, which get damped fast by 
direct dissipation (presumably through shock formation). 

Our calculations support the conjecture given by eq. Q. 
According to it the coupling of Alfven and compressible mo- 
tions is important only at the energy injection scales where 
5Vi ~ Va- As the turbulence evolves the perturbations be- 
come smaller and the coupling less efficient. Typically for 
numerical simulations the inertial range is rather small and 
this could explain why marginal coupling of modes was not 
noticed. 

Our results show that MHD turbulence damping does 
not depend on whether the fluid is compressible or not. The 
incompressible motions damp also within one eddy turnover 
time^°. This is the consequence of the fact that within the 
strong turbulence^^ mixing motions perpendicular to mag- 
netic field are hydrodynamic to high order (CLV02a) and 

with shocks, which are absent in incompressible simulations by 
Miiller & Biskamp (2000). 

It is generally believed that decay of turbulence energy follows 
a power-law. A possible expression is E{t) ~ Ea/[1 + C{t — fo)]" 
(C. McKee, private communication), where Eq is the energy at 
the initial time to. Our claim is that decay of turbulence follows 
E{t) ~ -Eo/[l-|-C'(t-io)/tcas,o]"> where tcas.o is the eddy turn- 
over time at t = to and C' is a constant of order unity. 

For a formal definition of strong, weak and intermediate tur- 
bulence, see Goldreich & Sridhar (1997) and CLV03a, but here 
we just mention in passing that in most astrophysically important 
cases the MHD turbulence is "strong" . 



the cascade of energy induced by those motions is similar to 
the hydrodynamic one, i.e. energy cascade happens within 
an eddy turnover time. 

The reported (see Mac Low et al. 1998) decay of the 
total energy of turbulent motions Etot follows t~^ which can 
be understood if we account for the fact that the energy 
is being injected at the scale smaller than the scale of the 
system. Therefore some energy originally diffuses to larger 
scales through the inverse cascade. Our calculations stim- 
ulated by illuminating discussions with Chris McKee show 
that if this energy transfer is artificially prevented by inject- 
ing the energy on the scale of the computational box, the 
scaling of Etot becomes closer to t~^ . 

Can compressible MHD turbulence decay 
slowly? Incompressible MHD computations (see Maron & 
Goldreich 2001; CLV02a) show that the rate of turbulence 
decay depends on the degree of turbulence imbalance^^, i.e. 
the difference in the energy fluxes moving in opposite direc- 
tions. The strongly imbalanced incompressible turbulence 
was shown to persist longer than its balanced counterpart. 
This enabled CLV02a to speculate that this may enable en- 
ergy transfer between clouds and may explain the observed 
turbulent linewidths of GMCs without evident star forma- 
tion. Imbalanced turbulence can also make it possible to 
transfer energy from the galactic disk to heat the Reynolds 
layer (see Reynolds 1988). 

Our results above show a marginal coupling of com- 
pressible and incompressible modes. This is suggestive that 
the results obtained in incompressible simulations are ap- 
plicable to compressible environments if amplitudes of per- 
turbations are not large. The complication arises from the 
existence of the parametric instability (Del Zanna, Velli, & 
Londrillo 2001) that happens as the density perturbations 
refiect Alfven waves and grow in amplitude. This instability 
eventually controls the degree of imbalance that is achiev- 
able. However, the growth rate of the instability is substan- 
tially slower than the Alfven wave oscillation rate. There- 
fore, if we take into account that interstellar sources are 
intermittent not only in space, but also in time, the trans- 
port of turbulent energy described in CLV02a seems feasible. 
Here we mention that the growth of the parametric instabil- 
ity described above may provide an alternative explanation 
for the observed infall motions of the cloud cores (Tafalla 
et al. 1998). Earlier those motions were explained as aris- 
ing from linear damping of Alfven waves (Myers & Lazarian 
1998). 

Is the correlation between magnetic field and 
density tight in MHD turbulence? In the traditional 
static paradigm of the ISM, density and magnetic field in- 
crease simultaneously as clouds contract. Introduction of 
turbulence in the picture of ISM complicates the analysis 
(see discussion in Vazquez-Semadeni et al. 2000; CLV03a). 
Our results confirm earlier claims (e.g. Passot & Vazquez- 
Semadeni 2003) that magnetic field - density correlations 
may be weak. First of all, some magnetic field fluctuations 
are related to Alfvenic turbulence which does not compress 
the medium. Second, slow modes in low /3 plasmas are es- 
sentially density perturbations that propagate along mag- 



This quantity is also called cross helicity (see Matthaeus, Gold- 
stein, & Montgomery 1983). 
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Figure 10. /3 = 1 case. {a)left: Alfvcn mode spectra, {b) second-left: Slow mode spectra, {c) second-right: Fast mode spectra. {d)right: 
Density structure. 



netic field and which marginally perturb magnetic fields. 
Moreover, our calculations (see Fig. |5J| show that at sub- 
stantially high Mach numbers the density correlations do 
not show anisotropics, while magnetic field fluctuations are 
anisotropic. On the basis of our calculations we might ex- 
pect that the correlation may be a bit higher for high /3 than 
for low /3 plasmas. 

Can viscously damped turbulence explain 
TSAS? The term "tiny-scale atomic structures" or TSAS 
was introduced by Heiles (1997) to describe the mysteri- 
ous H I absorbing structures on scales from thousands to 
tens of AU, discovered by Dieter, Welch & Romney (1976). 
Analogs of TSAS are observed in Nal and Call (Meyer & 
Blades 1996; Faison & Goss 2001; Andrews, Meyer, & Lau- 
roesch 2001) and in molecular gas (Marscher, Moore, & Ba- 
nia 1993). 

Those structures can be produced by turbulence with a 
spectrum substantially more shallow than the Kolmogorov 
one, e.g. with the spectrum E{k) ~ (see Deshpande 
2000). Simulations in CLV02b and theoretical calculations 
in Lazarian, Vishniac & Cho (2003) show that the magnetic 
field in the viscosity-damped regime of MHD turbulence (see 
|SJl can produce such a spectrum of magnetic fluctuations. 
Our calculations above are indicative that this will translate 
in the corresponding shallow spectrum of density. Our cal- 
culations are applicable on scales from the viscous damping 
scale (determined by equating the energy transfer rate with 
the viscous damping rate; ~ 0.1 pc in the Warm Neutral 
Medium with n = 0.4 cm"^, T= 6000 K) to the ion-neutral 
decoupling scale (the scale at which viscous drag on ions 
becomes comparable to the neutral drag; ^0.1 pc). Be- 
low the viscous scale the fluctuations of magnetic field obey 
the damped regime shown in Figure|5j3 and produce density 
fluctuations. For typical Cold Neutral Medium gas, the scale 
of neutral-ion decoupling decreases to ~ 70 AU, and is less 
for denser gas. TSAS may be created by strongly nonlinear 
MHD turbulence! 



6.3 Application of the scaling laws obtained 

Many astrophysical problems require some knowledge of the 
scaling properties of turbulence. Therefore we expect a wide 
range of applications of the established scaling relations. 
Here we discuss how our understanding of MHD turbulence 
affects a few selected issues. 

As we mentioned in the introduction, the observations 
(see CLV03a for a review) indicate that interstellar spec- 
trum exhibits Kolmogorov-type scaling E{k) ~ k~^^^ . On 
the basis of what we learned by now (see also Higdon 1983, 



GS95, Lithwick & Goldreich 2001, CL02) we may conclude 
that the MHD spectrum is different from the spectrum 
of isotropic Kolmogorov turbulence. Anisotropy of Alfvenic 
and slow modes ensures that the observationally measured 
spectrum is dominated by perturbations perpendicular to 
the magnetic field, i.e. E(k) ~ k'^^^ . Admixture of fast 
modes can alter the spectrum, but observations indicate that 
those modes do not dominate the detected signal. 

Whether or not one can use Kolmogorov isotropic scal- 
ing for practical calculations depends on the nature of the 
problem. For instance, Cho & Lazarian (2002b) use Kol- 
mogorov scaling to calculate the spectra of fluctuations aris- 
ing from synchrotron emission and of fluctuations of the de- 
gree of starlight polarization and obtained good correspon- 
dence with the observed statistics. However, for many prob- 
lems turbulence anisotropy is essential. 

Cosmic rays. — The propagation of cosmic rays is 
mainly determined by their interactions with electromag- 
netic fluctuations in the interstellar medium. For practi- 
cal calculations it is usually assumed that the turbulence 
is isotropic with a Kolmogorov spectrum (e.g., Schlickeiser 
& Miller 1998). How should these calculations be modified 
in view of our findings? Chandran (2000) and Yan & Lazar- 
ian (2002; see also Lazarian, Cho & Yan 2002 for a review) 
calculated the efficiency of cosmic ray scattering by Alfvenic 
modes and found that, if the energy injected at large scales, 
the scattering by Alfvenic modes is absolutely negligible. 
The difference between the calculations in Yan & Lazarian 
(2002), which used the tensor description of Alfven waves ob- 
tained in CLV02a, and the calculations that use Kolmogorov 
turbulence was 10 orders of magnitude! Yan & Lazarian 
(2002; 2003a) identified fast modes as the principal source 
of cosmic ray scattering and used the CL02 results to calcu- 
late the scattering rate in low /3 plasma. Damping of those 
modes, however, depends on the value of /3, which entails 
the dependence of cosmic ray scattering on the temperature 
and density of the background plasmas. 

Dust dynamics. — Turbulence induces relative dust 
grain motions and leads to grain-grain collisions. These col- 
lisions determine grain size distribution, which affects most 
dust properties, including absorption and H2 formation. Un- 
fortunately, as in the case of cosmic rays, earlier work ap- 
pealed to hydrodynamic turbulence to predict grain rela- 
tive velocities. Lazarian & Yan (2002) and Yan & Lazarian 
(2003b) considered motions of charged grains in anisotropic 
MHD turbulence and found that the direct interaction of 
the charged grains with turbulent magnetic field results in 
a stochastic acceleration that can potentially provide grains 
with supersonic velocities. 
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Other applications. — The obtained scahng laws are 
essential for understanding the density fluctuations within 
HII regions (Lithwick & Goldreich 2001), stochastic mag- 
netic reconnection in fully ionized (Lazarian & Vishniac 
1999, 2000) and partially ionized plasma (Lazarian, Vish- 
niac & Cho 2003), for energy transfer to electrons in gamma 
ray bursts and solar flares (see Lazarian, Petrosian, Yan, & 
Cho 2003). This list can be easily extended. 



7 SUMMARY 

In the paper, we have studied statistics of compressible MHD 
turbulence for high, intermediate, and low f5 plasmas and for 
different sonic and Alfven Mach numbers. For subAlfvenic 
turbulence we provided the decomposition of turbulence into 
Alfven, slow and fast modes. We have found that the cou- 
pling of compressible and incompressible modes is weak and, 
contrary to the common belief, the drain of energy from 
Alfven to compressible modes is marginal along the cascade. 
For the cases studied, we have found that GS95 scaling is 
valid for Alfven modes: 

Alfven; E^[k) oc oc fc^/^ 

Slow modes also follows the GS95 model for both high /3 and 
mildly supersonic low /3 cases: 

Slow: E'ik) (xk-^^^, k^\ (X k]^'-' . 

For the highly supersonic low /3 case, the kinetic energy spec- 
trum of slow modes tends to be steeper, which may be re- 
lated to the formation of shocks. 

Fast mode spectra are compatible with acoustic turbulence 
scaling relations: 

Fast: (k) oc k~^^^ , isotropic spectrum. 

Our super- Alfvenic turbulence simulations suggest that the 
picture above holds true at sufficiently small scales at which 
Alfven speed Va is larger than the turbulent velocity vi . This 
part of our study shows that compressible MHD turbulence 
is not a mess. On the contrary, its statistics obeys well deter- 
mined universal scaling relations. The importance of these 
relations for different branches of Astrophysics is obvious. 

Addressing the issue of MHD turbulence damping in 
partially ionized gas we showed that the viscosity-damped 
regime of MHD turbulence below the viscous damping scale 
that was reported in CLV02b for incompressible fluid per- 
sists when compressibility is present. The spectrum of den- 
sity fluctuations that we obtain may be related to the mys- 
terious tiny scale structures observed in the ISM. 
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APPENDIX A: MODE DECOMPOSITION 

Let us consider a small perturbation in the presence of a 
strong mean magnetic field. We write density, velocity, pres- 
sure, and magnetic field as the sum of constant and fluc- 
tuating parts: p ^ po + p, v ^ vq + v, P ^ Pq + p, and 
B —> Bo -I- b, respectively. We assume that vo = and that 
perturbation is small : p « po, etc. Ignoring the second 
and higher order contributions, we can rewrite the MHD 
equations as follows: 

|^+poV-v = 0, (Al) 

po^ + V{a^p) - -L(V X b) X Bo = 0, (A2) 
at Atv 

5 + V X [v X Bo] = 0, (A3) 
ot 

where we assume a polytropic equation of state: p — p 
with = 7po/po- We follow arguments in Thompson (1962) 
to derive magnetosonic waves. Let ^(r, t) be the displace- 
ment vector, so that d^/dt — v. Assuming that the dis- 
placements vanish at t = 0, we can integrate the equations 



as follows 

p + poV ■ C = 0, (A4) 

C = aV(V-C) + (V X b) X Bo/47rpo, (A5) 

b = V X (5 X Bo). (A6) 



The momentum equation (eq. IA5t becomes 

i = aV(V ■ C) + [V X (V X (e X Bo))] X Bo/47rpo 
= aV(V ■ + V(Bo V • C - Bo ■ VBo • 0/4npo 

- (Bo ■ Vf^/4npo + [Bo(Bo ■ V)V ■ f] /47rpo (A7) 

Using a = o^/Va = /9(7/2), Va = Bo/47rpo, we have 

C/Kl-V[(a4-l)V-C-(Bo-V)(Bo-C)]-(Bo-V)'C+(Bo-V)(V-OBo = 0(A8) 

In Fourier spacethe equation becomes 

i/v2 + fck[(Q + l)fcCfc - + ~ hKkh = 0, (A9) 

where Cfc — ? ' k, = ^ • ky, k = k/fc, and ky is unit vector 
parallel to Bo (i.e. ky = Bo). Assuming^ = —ui'^S^ = —c^k^S^, 
we can rewrite IIA9t as 

(cVv|-cos^6i)5-[(Q-Hl)5fe-cos6'C||]k-fcos6i5fefc|l = 0, (AlO) 

where cos6' = k\\/k and is the angle between k and Bo. 
Using k = sin Oli_i_ + cos Sky , we get 

(cVKI - cos^ 6*)^ - [(a -f- l)ffe - cose^y] sin6'kx 
-{[(a + l)Cfc-cos6l^y]cos6'-cos6'^fc}ky =0. (All) 
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Writing ^ = ^±k_L + + Cv<^i we get 

[c/Vl ~ cos^ 6l)5x - [(a + 1)6 - cos e*^!,] sin^ = 0, (A12) 

{c/Vl - cos^ 6i)C|| - Kfc - cose^ii] cosf = 0, (A13) 

(c/Vl - cos" e)^^ = 0. (A14) 

The non-trivial solution of equation IIA14t is the Alfven 
wave, whose dispersion relation is Lj/k — VacosS. The di- 
rection of the displacement vector for Alfven wave is parallel 
to the azimuthal basis ip: 

|a = -(^ = kx X k||. (A15) 

Let us consider solutions of equations HA 1211 and 1A13II . 
Using 6 = 5x sin 6 + cos 6, we get 

(cVVa - cos^ - (q + 1) sin^ — Q COS 8 sin 6'^|| — 0j(A16) 
{c^/Va - cos^6l)5|| - Q sin e cos e^x - (q - 1) cos^ 6IC|| = 0(A17) 
Rearranging these, we get 

{c^/Va - Qsin^e - l)^x - Q cose sin 1 1 = 0, (A18) 
{c^/Va - a cos^ 61)^11 - asinecose^x = 0. (A19) 
Combining these two, we get 



The slow basis ^ lies between ky and ~9. The slow basis 
lies between kx and k (Fig. 0. Here overall sign of ^ and 
is not important. 
When a ^ 0, equations 1A30II and llA29l l becomes 



6 ~ k|| — (q sin 5 cos 6')kx, 
(,f ~ (q sin 61 cos 6')k|| -I- kx. 



(A31) 
(A32) 



In this limit, ^ is mostly proportional to ky and ^/ to ki 
When Q — > oo, equations <A30t and 1A29II becomes 



6 ^ -61+ (sin6'cose'/Q)k, 
if ^ {sinecos9/a)e + 



(A33) 
(A34) 



(c /Va — a sin 6 — l){c /Va — a cos 

2 ■ 2 o 2 

= a sm a cos 
Therefore, the dispersion relation is 
c'^/Vi - (1 + a)cVvi acos^ 6^0. 
The roots of the equation are 

(1 + a 



± ^(l + a)2-4ac 



(A20) 
(A21) 

(A22) 



2 1 , r2 

where subscripts 'f and 's' stand for 'fast' and 'slow' waves, 
respectively. 

We can write 



^ = ^||k|| -f ^kx oc 



C||fcx 
Cxfcii 



fciikii + fc I k I 



(A23) 



Plugging eq. 1A22II into eq. HAlSt and (IA19II . we get 



1 + Of , 



2 2 

1+a 



± a cos 

2 



where D — {1 + a)^ 
k± = k cos 9, we get 



Cx = Q cos 6* sin 6*^11, (A24) 

^11 = a cos 6* sin 6'^x , ( A25) 
4acos^&. Using kn — kcos9 and 



-1 



2 

1+a 



2 2 
Arranging these, we get 
Cllfcx -l + a±VO 



^xfc|| — asin^S^xfc|| = acos^ 6'^||fcxj(A26) 



C||fcx -Qcos^ e^iifcx =asin^ 6lCxfc||(A27) 



(A28) 



where the upper signs are for fast mode and the lower signs 
for slow mode. Therefore, we get 

^ ^ (A29) 

(A30) 



6 oc (-1 + a - V-D)fc||k|| + (1 + Q - 
if oc (-1 + Q + \/D)fc||k|| + (1 + a - 



D)fcxkx, 
"D)fcxkx. 



When Q = oo, slow modes are called psewiio- Alfvenic modes. 

We can obtain slow and fast velocity component by 
projecting Fourier velocity component Vk onto ^ and 
respectively. 

To separate slow and fast magnetic modes, we assume 
the linearized continuity equation {ijpk = pok ■ vj,) and the 
induction equation {uihk = k x (Bo x vj,)) are statistically 
true. From these, we get Fourier components of density and 
non-Alfvemc magnetic field: 

pk = (poAufc,s/cs)k ■ 6 + (/9oA?;fej/c/)k • 6 

= Pk,s + Pk,f, (A35) 

bk = {BoAvk,s/cs)\Bo X is\ + {BoAvk,f/cf)\Bo X if\ 

= bk,s + bk,f, (A36) 

= Pk,s{Bn/pn){\Bo x ^^l/k ■ is) 

+ pfc,/(Bo/po)(|Bo xOI/k-C/), (A37) 

where Avk (x — (superscripts '+' and '-' represent 
opposite directions of wave propagation) and subscripts 's' 
and 'f stand for 'slow' and 'fast' modes, respectively. From 
equations 1A35II . 1A36II . and llA37ll . we can obtain pk^s, pkj, 
bk,s, and bkj in Fourier space. 



APPENDIX B: SCALE-DEPENDENT 
ANISOTROPY 

Fig. lBlh shows the shapes of Alfven eddies of different sizes. 
Left 3 panels show an increased anisotropy as we move from 
the top (large eddies) to the bottom (small eddies) . The hor- 
izontal axes of the left panels are parallel to Bo. Structures 
in the perpendicular plane (right panels) do not show a sys- 
tematic elongation. However, Fig. IBlb shows that velocity 
of fast modes exhibit isotropy. Data are from a simulation 
with 216^ grid points, Ms = 2.3, and (3 = 0.2. 



© 0000 RAS, MNRAS 000, 000-000 



compressible MHD 17 




Figure Bl. Anisotropy as a function of scale. {Left) Alfven mode velocity show scale-dependent anisotropy. {Right) Fast mode velocity 
show isotropy. Only part of the data cube is shown. Lighter tones are for larger |v|. Magnetic field show similar behaviors. 
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